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Abstract 



We present a strategy for solving time-dependent problems on grids with local re- 
finements in time using different time steps in different regions of space. We discuss 
and analyze two conservative approximations based on finite volume with piecewise 
constant projections and domain decomposition techniques. Next we present an itera- 
tive method for solving the composite-grid system that reduces to solution of standard 
problems with standard time stepping on the coarse and fine grids. At every step of the 
algorithm, conservativity is ensured. Finally, numerical results illustrate the accuracy 
of the proposed methods. 
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1 Introduction 



In many physical applications, there are special features which greatly affect the solution 
globally as well as locally. One important example is the local spatial and temporal behavior 
of multiphase fluid flow around a production well in the petroleum recovery applications. 
To capture this local behavior, spatial local refinement is necessary. However, it requires a 
reduction of the time step, compared to the one used with a coarse mesh, in order to get 
a solution accurate enough in the refined zone and to avoid convergence problems when 
solving the non linear discretized equations. When applied uniformly on all the simulation 
domain, this reduced time step leads to unacceptable cpu-time making the use of local time 
steps highly desirable. To be efficient, a local time-stepping strategy (numerical scheme 
and solution method) must 

• ensure accuracy of the solution i.e. the solution has to be more accurate than the 
one obtained with a global coarse mesh, 

• ensure stability without any too restrictive condition on the time step, 

• lead to reduced cpu-time compared to the one obtained when using a small time step 
on the whole domain. 

In the framework of reservoir simulation where local grid refinement is necessary to repre- 
sent correctly important local phenomena in the wells vicinity, the corresponding numerical 
scheme must also be locally conservative in order to be applicable to multi-phase flow sim- 
ulations where a coupled system of parabolic and hyperbolic equations has to be solved. 

For parabolic equations, different approaches have been proposed in the past which 
extend the classical implicit finite difference scheme to local refinement in time. In [6], 
the scheme is written as a cell centered Finite Volume scheme. At the interface between 
the coarse time step zone and the refined time step one, the flux over the coarse step 
is taken equal to the integral over the corresponding refined steps of the flux computed 
from the refined zone. This refined flux approximation requires values of the unknowns 
in the coarse zone at small time steps which are computed using piecewise constant or 
linear interpolation from the coarse unknowns. The scheme is conservative. Stability and 
error estimation are obtained for the piecewise constant interpolation. On the contrary, 
for the linear interpolation, stability is obtained under a sufficient condition which is as 
restrictive as the time step limitation obtained for an explicit scheme. In [3], Dawson et 
al proposed to couple classical implicit finite difference schemes in the refined and coarse 
zones using an explicit approximation at the interface on a larger mesh size in order to 
attenuate the time step limitation due to the explicit approximation. Although interesting 
for its simplicity, this approach can not be retained due to its time step limitation. In [5], 
Ewing and Lazarov proposed an implicit non conservative approach. The scheme for the 
coarse nodes is straightforward while the fine grid nodes located at the interface between 
the coarse and fine regions require "slave" points at small time steps on the coarse grid 
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side, which are not grid points. As in [6], the values of the unknown at these slave points 
are obtained by linear interpolation in time between the corresponding nodes of the coarse 
grid, and the set of discretized equations involves all the unknowns between two coarse time 
levels. Stability and error analysis are performed. The solution method uses an iterative 
method associated to a coarse grid preconditionner of the Schur complement of the system 
where the refined region unknowns have been eliminated. In the more applied framework of 
compositional multiphase flow, [3] introduced an implicit timestepping method. For each 
global time step, the problem is solved implicitly in the whole domain but using a linear 
approximation of the model in the refined regions which avoids any convergence problem 
of the non linear solver due to refined mesh. Then, the refined zones are solved using a 
local time step and taking as boundary conditions the fluxes computed during the first 
stage at the interface between the refined and coarse zones. This approach ensures that 
the method is conservative. It is moreover rather efficient as, compared to the cpu-time 
necessary to solve the problem with a large timestep on the whole domain, it only requires 
additional cpu-time to solve the equations once in the refined zone. However, the accuracy 
of the solution is not controlled. Looking for an efficient solution method, [12\ used the 
same finite difference scheme as [5j for linear parabolic equations but proposed a predictor- 
corrector method. In the predictor stage, the solution is computed at the coarse time step 
on all the domain and in the correction step, the solution is computed in the refined grid 
at small time steps using values at slave nodes interpolated from the coarse nodes solution 
obtained in the first stage. They show that the predictor corrector approach preserves the 
maximum principle satisfied by the solution of the scheme. 

Our paper proposes a local time step strategy based on the domain decomposition 
framework. It extends the approach introduced in [6] by generalizing the interface condi- 
tions used to couple the coarse and refined time-step domains. The method is conservative. 
Stability and error estimates, which are different from that obtained in [6], are presented. 
A solution method which improves the predictor-corrector methods of [1] and [12] is pro- 
posed. In order to simplify the presentation and to concentrate on the difficulties arising 
from the local refinement in the time direction, we will first explain the approach in the 
case of a one-dimensional spatial problem in section [2j Stability and error estimates are 
then proven in the more general case of nD spatial grids in sections 3 to 6. Finally, some 
numerical results are presented in section 7. 

2 Description of the local time stepping strategy 

We consider the following problem: Let T > and f2 be an open bounded domain of M. d , 
d ^ 1, po : O i->- R and / : O X (0, T) i-> R be given functions. Find p : Q x [0, T) i-> R such 
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that: 

^r(x,t) - Ap(x,t) = f(x,t) VxG0 Vte[0,T] 
at 

p(0, x) = po(x) Vi G 
p(x, t) = Vx G 80 Vi G [0, T] 

In order to explain the scheme, we consider the d = 1 case with the cell centered grid 
shown in Figure Q] and a time step which is variable in space. Namely, the domain is 
decomposed into two non overlapping subdomains f2i and VL2 where two different time- 
step sizes are used : the coarser time step is denoted 5t2 (in O2) and the finer time step is 
denoted 5h (in ft x ) such that fZSh = 5t 2 with JC G N* (Figure 
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2.1 Discretization 



In each sub domain, the equation is discretised using a classical cell centered finite volume 
implicit scheme: 



|L( p ^+l _ p f) _ (tt £H-l _ ^+1) = ^ jn.H-1 for < j (2a) 

^r 1 - - = for * > j ( 2b ) 

where is an approximation of the unknown in the space-time cell (x ■ 1 , x • , 1 ) x (t 1 , t , 1), 
p"' in the cell (aJ J _i , i ) x (t n fe _i , t n fc+ i ) (see Figures [Q and ED , 



' 3 /ij-eJti 
and 



x . i ./< i 
2-3 ™-? 



U = TTT / / f(x,t)dxdt. 
■> hjdt 2 J x . 1 J t 1 

Except for the boundary nodes (which are handled classically and not precised here) and 
interface node x I+ i, the flux approximation is given by: 

p^. 

u s . , 1 = — for s = n if 7 > J or s = n, A; if 7 < J (3) 

For the approximation of the fluxes on the interface and it?' 1 , we consider the 

space-time domain decomposition framework and introduce p". + \ and p™' k i the unknown 

+ 2 + 2 

approximations on the interface x I+ i. The flux approximations are classically obtained as 

n.k n,k 
P I+ l~Pl 

n,k > 2 I a\ 

u I+ i = — ( 4 ) 
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0.5/i/ 



n+l n+l 

u n+l = 1+2 (5) 

I+h 0.5/iz+i v ; 

Discretizations in the two domains are linked by interface conditions which enforce flux 
and unknown continuity on the interface x I+ i x (t n _i,t n+ 1). We consider the following 
two possibilities: 
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K 



5 ^ u lX\ = Y, 5t ^ n iii ( 6a ) 

fc=l 



P n T '\=p n T tl fc = l--JC (6b) 



or: 



'+2 ^+2 



^; fc i=<i fc = l---/C (7a) 

2 2 

These interface conditions can be rewritten in terms of the L2 orthogonal projections on 
sets of piecewise constant functions in time (see section I3.3[) . Both sets of conditions ensure 
local conservation for the coarse time step. 

Another way to couple the fine and coarse grid which is somewhat more natural, is to 
introduce the unknown approximations p™ + \ and not on the interface but rather in 

the neighbouring cells p™ +1 and Pj+i- The flux approximations on the interface are then 
directly expressed as for an interior edge and instead of (HD,©, we have: 

n,k n,k 

n,k Pl+1 ~~ Pi , s 

U M = h I+1 (8) 
n+1 _ n+1 

«£1 = Pl+ \ Pl (9) 

h I+h 

and the interface conditions become: 

K 



or: 



2 i% 2 

pT+!=Pltl k = !■■■)€ (10b) 



«?+i= u £| k = l--K (11a) 

st 2P ^ +i =j2^p? k ( llb ) 

k=l 



We have thus four possible coupling schemes: . - — - - <|7|> — <|B|> - - <fT0|> — <|g|> . <|g|l - <PCT] 
which are analysed in the sequel. Equations ([8]) . ([9]) . fp~Q|) are the approximations proposed 
in [6]. 
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2.2 Solution method 



To solve the system of algebraic equations for the unknowns values of the approximate 
solution between two coarse time levels, which includes all intermediate local time levels, we 
propose a method which combines the attractive feature of predictor-corrector approaches 
with the accuracy of domain decomposition type iterative algorithms. The method includes 
a predictor stage which corresponds to the computation of the solution on the coarse 
grid time and an iterative corrector stage where refined and coarse grids unknowns are 
solved alternatively, until interface conditions are satisfied, using a Schwarz multiplicative 
Dirichlet Neumann algorithm [11]. If we consider (f8]).([9j). (|10p interface conditions, the 
algorithm is: 

• Predictor stage: computation of an approximate solution at coarse time step on the 
whole grid: P™ +1 for all j 

^(p7 +1 - p?) - - = vr 1 fOT a11 j ( 12 ) 

• Corrector iterative stage if interface conditions (jlOp are used: 

Solve alternatively the equations in domain Q\ using (llObl) interface condition and in 
domain using (jlOap . until both interface conditions are satisfied simultaneously. 

If ([8j).(i9l). (fTTj) interface conditions are used, the corrector iterative stage is : 

• Corrector iterative stage: 

Solve alternatively the equations in domain 17 1 using (|llap interface condition and in 
domain using (lllbp . until both interface conditions are satisfied simultaneously. 

These algorithms can also be written for (HI),©,© or ([I]),©,© interface conditions. We 
can notice that it is not necessary to iterate the corrector stage until convergence to obtain 
a conservative solution. It is sufficient to stop the process after the resolution of the 
domain where Neumann interface conditions are imposed. The algorithm proposed in [4] 
consists in the predictor stage (|12p and in the first iteration of the corrective stage with 
(lllap interface condition. This solution method is not limited to linear discrete equations 
and can be extended to the resolution of the non linear equations that arise in petroleum 
recovery applications. Following the idea of [4j, the predictor stage would then use a linear 
approximation of the problem in the refined domain in order to avoid convergence problems 
of the Newton algorithm, while the iterative corrector stage would consider the non linear 
problem. 



S 



3 Finite volume discretization 



Problem (Q} is rewritten as a domain decomposition problem. The domain 0, is decomposed 
into two non overlapping subdomains Qi and O2 U O2 = H and Q\ fl O2 = 0)- The 
interface is denoted by T = Oi n O2. Problem (P) is equivalent to: 
Find p 1 : Oi x [0, T] h-> R and p 2 : tt 2 x [0, T] i-> R such that: 

^-(a?,t)-Ap i (a;,t) = /(aj,t) VxGO, Vie[0,T] Vie {1,2} (13a) 

p*(z,0) =po(a:) Vx G fi t Vi€{l,2} (13b) 

p i (x,t) = o yxedttindn v*g[o,t] v«e {1,2} (13c) 

^-(x,*) + ^-(a?,t) = VxGT Vte[0,rl (13d) 
p 1 (x,t) =p 2 (x,t) VxeT Vie [0,2*] (13e) 

where raj is the outward normal to domain f2j, z = 1,2. 

Problem ()13p is discretized using a cell centered finite volume scheme in each subdomain 
[7]. We choose this scheme as an example but other schemes would be possible as well. 

3.1 Mesh and definition 

For i = 1,2, let % be a set of closed polygonal subsets associated with such that 
Cti = LixeTiK. We shall denote h = max^e{i J 2},K'e7j diam{K) its mesh size. We shall use 
the following notation for all i = 1,2. 

• £q { is the set of faces of %. 

• £io is the set of faces such that <9f2j fl Oil = {J el =£ iD e (let us recall that a Dirichlet 
boundary condition will be imposed on dtti n dQ). 

• £i is the set of faces such that dVti\d£l = Ujg^e. Grids are matching on the interface 
so that 

S:=£ 1 = £ 2 . 

• vk e%, 

£{K) denotes the set of faces of K. 

£id{K) = £(K) n £iD is the set of faces of K which are on d£li n d£l. 

£i{K) = £{K) n £i is the set of faces of K which are on d£li\dVt. 

Ni{K) = {K 1 e %\K (1 K' € £n t } is the set of the control cells adjacent to K in 

Ki(e) denotes the cell of % adjacent to e € £{ U £j£>. 



9 



Figure 3: Assumption 13.11 



• the time step in subdomain 0j is denoted by 6ti, and N% denotes the number of time 
steps of the simulation so that N\5ti = = T. Parameters 5U, Ni satisfy 

• Let [0, T]sn denote the discretization of the time interval [0, T] in each subdomain 
Ui : [0,T]st t = (t l n )n=i,...,Ni, with t % n = (n — \)5ti ; since the time discretization in 
ill is a refinement of that in we shall also write : k = (n — 1)^2 + (k — ^)St\, 
n = 1, ... ,N 2 , k = 1, ... ,/C 

We make the following geometrical assumptions on the global mesh : T = T\ U T 2 

Assumption 3.1. T is a finite volume admissible mesh, i.e., T is a set of closed subsets 
of dimension d such that 

• for any (K, K') G T 2 with K / AT', one has either [K K'\ := K n K' G % U Sq 2 or 
dim(AT n if') < d - 1 

• for i = 1,2, there exist points (y e ) eg £ f2 on the faces and points {xk)k^% inside the 
cells such that (see figure [3]) 

— for any adjacent cells K and K', the straight line [xk,xk>] is perpendicular to 
the face [K K'\ and [x K ,x K >\ n [AT K'] = {y e } 

— for any face e G £j£>, let K(e) € % be such that e C AT: then the straight line 
[ x K(e) j is perpendicular to e 

• Each mesh i = 1,2 has at least one interior cell. 



Notation 3.2. For all e G Si UShj, i = 1, 2, denotes the distance between x^^e) an d y, 
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3.2 Cell centered finite volume scheme in the subdomains 

The unknowns of the scheme and what they aim to approximate are (i = 1, 2) 

v 1 k - p l (x K ,t n ), K e % 

u l t ~ ^-{ye,t n ), e G £iU£ iD 
The numerical flux is defined by: 

i,n i,n 



u % f = l{£) Ve E £ U £ iD Vn G {O..Ni - 1} Vi G {1, 2} (15) 



The scheme is defined by (see e.g. [7] for its derivation): 

i,n+l i,n i,n+l i,n+l 

Pk ~PK meas{K) _ y Vk> _Pk meas([ ^]) 
SU K J ^ d(x K ,x K i) u Jy 

^ ?4' n+1 meas(£) = /^ n+1 meas(K) V-ftT G % Vn G O.JV; - 1 (16) 

ee£i(K)L)£ iD (K) 

where cI(xk,xk') is the distance between xk and xk 1 and /,£ is an approximation to 
3F J^ n+1/2 meas(_ft:) /if /• ^he initial and boundary conditions are discretized by: 



- 

J,0 



pk=Po(a*:) VKG^ Vi€{l,2} (17) 
p i,n+i = o Ve G «S, D Vn G 0. JVj - 1 G {1, 2} (18) 

When there is no domain decomposition, this scheme has been analyzed in [7] in the more 
general case of discontinuous coefficients, and it is proven to be of order 1 for a discrete 
i? 1 -norm. 

In order to define the domain decomposition discretization scheme, we shall define in 
section [33] the matching conditions for the diffusive fluxes. 
Discrete spaces 
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• Similarly, Pq{£ x [0, T]stJ is the space of piecewise constant functions on the interface 
for the time mesh of subdomain Qj. 

. Po([0,T] dti ) = \p:[0,T]^R\ Vn e {0, . . . , N { - 1} t i } = C* } 

^ lv n — 1/2' ra+l/2 y J 

These are spaces of piecewise constant functions. 

Let p l G Po(1i x [OjT]^), we denote its restriction 

• to Qi x {4} by ^' n , for n € {0, . . . , JVi - 1} 

• to e x [0, T] 5ti by for e G £ 

• to S X [0, T]^ by pi 

We introduce the following norms and semi-norms for any p l G Po(% x [0, T]^. ): 

ll^liw = E (^) 2 meas(^) 

Ni 

11^11^(0,^(00) = Ew 2 



and 



and 



i 2 (0i) 

ra=0 



(/ i,n i,n\2 
E T meas([^]) 

i i,n i,n\2 i i,n\2 

+ E M ~ PK meas ( £ )+ E #^meas( e ) 



(20) 



bii^^E^i^H^ ( 21 ) 

n=0 

Definition 3.1. Letp l ,u l G Po(^i x [0, Tl^J, * = 1>2, we define a discrete scalar product 
by: 

2 2 Ni-1 

E< u< 'P%(o,T ; L 2 (r)) =E E ^E M e" +1 Pr +1 ^a S (e) 

i=l i=l n=0 e££ 

Notation 3.3. Let i = 1,2, p l G Po(%, X [O,? 1 ]^.), /or e G f , ?4(p) denotes the associated 
numerical flux defined by ()15p . Very often, we will simply write u\ and Ug = (u e (p l )) ee £ i . 
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3.3 Finite volume on the interface 



In order to enforce the weak continuity of the primary unknown p and of its normal 
derivative (denoted by u) across the interface F x [0, T], we introduce Qi the L 2 orthogonal 
projection onto Po([0, T]^). We have the following compatibility condition: 

Lemma 3.4. For all m € P ([0,T] sti ), i = 1,2, 

(Qi(u2),Ui)l2([0,T]) = {U2, Q2(«l)>L2([0,T]) 

As in mortar methods [2], we consider that one subdomain enforces the weak continuity 
of the primary unknown which is interpreted as the Dirichlet interface condition. This 
subdomain is called the master. The other subdomain is called the slave. It enforces 
the weak continuity of the normal derivative which corresponds to a Neumann interface 
condition. 

Since here the interfaces are non matching only in the time direction, it is possible to 
define matching conditions locally on each interface face e € £. 



Interface Scheme (ISi ) based on interface unknowns The subscript m will denote 



the master subdomain and e the slave ({m, e} = {1, 2}), the interface conditions on Tx [0, T] 

1 Pi =*«•) V£ ' ' (ISl) 

Overlapping interface scheme QIS2D The Dirichlet boundary condition is modified 
but not the Neumann one: 



p e e + d?ul=Q e (p?-d?u?) 



\/ee£ (IS 2 



The modified Dirichlet interface condition comes from the following relations: 



(pn _|_ 

Q e {p™ £ -d™ £ uf)-pt + dtu 



df + d% 

where d\ = d{xj^.^,y £ ). The first line of the above formula is somewhat natural. When 
writing the finite volume scheme for a cell K e (e) adjacent to the interface in the "slave" 
subdomain, it is necessary to approximate the flux on the face e. This is done using 
pressure values from both sides of the interface: the pressure in the "slave" subdomain 
and pressures values in the neighboring "master" subdomain. These last values are made 
compatible with the "slave" unknowns by using the transmission operator Q e . Finally, all 
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quantities are expressed in terms of interface values in order to ease a comparison with 

(USD. 

Due to the fact that the large time step 5t 2 is a multiple of the small time step 5tx, we 
have a simple form for the L 2 projection operators. 

Lemma 3.5. We have Q 2 : P o ([0,T] Stl ) ^ P ([0,T} St2 ) and Qx : P {[0,T} St2 ) h- P o ([0,T] Stl 
{ Forv 2 eP ([0,T} St2 ), Qi(v 2 ) m t i i/2) =^ 2 | (tl i/2 , t i i/a) VnGjO,...,^-!} 

I ^ n— 1/2' 71 + 1/2'' IV n— 1/2' n+1/2-' 

ForvtEPodO,^), «2 2 («i)) Kt 2 2 ) = -L /' n+1/ % 1 VnG{0,...,iV 2 -l} 

IV n-1/2' n+1/2' /, 2 

(22) 

We also need the following technical assumptions. For the family of meshes we consider, 
the mesh close to the interface is not too stretched: 

Assumption 3.6. There exists a constant a > such that d™ ^ ad e e for all e G £ ■ 
We also need a geometric assumption: 

Assumption 3.7. For all e G £, y £ is the barycenter of the face e and for i = 1,2, 

diamine) 2 



d(x Ki{£ ),y £ 

4 A general stability result 



0(h) 



We prove a stability result for (|16p modified by the introduction of discretization error 
terms R^, R % kk> anc ^ as well as -^k™ which will be defined precisely in the sequel 
(ie{l,2»: 

(i,n+l _ i,n \ / i,n+l _ i,n+l \ 

- *H — w - ^ (%^r - — ( w> 

«' n+1 -^ n+1 )meas(e) = F^ n+1 meas(if ) Vif G % Vn G {0, . . . , ty-l} . 

E e£i(isr)uf i£) (X) 

(23) 

Formula (|15p is modified as well by introducing error terms at each time step n 

i,n i.n 

U i,n = u i,n_ e Kde) y £e£ . u£w V nG{l,...,^} Vi € {1, 2} (24) 

a £ 

(with pl' n = for all s G fjo). 
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The interface conditions fllSiD are modified in the following manner by the terms V £ G 
Pq(£ x [0, T]$ te ) and U £ G Pq(£ X [0,T]^ m ) which will be defined in the sequel: 



u?=Q m {-u e e )+U? " 



Interface conditions (IIS2D are similarly modified for all s G £ : 



, Qe(PS m(e ))+^-^.( e) 
df + d% 

_ Q e { P T - d™uT + d™ur) + V £ -p e £ + d%u% - d%U t 



d™ + d e F 



The interface conditions ( IIS2D are thus modified in the following manner: 
( <=Q m (-0+Z4 



We make some additional assumptions. 

Assumption 4.1. i?^/ + R k >k = VK e% VST' G Ni{K) Vie {1,2} 

Assumption 4.2. VFe suppose that the following estimates are satisfied by the consistency 
errors: 

for all i G {1, 2} and all n G {1, . . . , Ni}, we have: 





=0{h^ 




\/k g % 


Vee£(tf)U 


K K 


=0(h^ 


+ K) 


MK G % 




ni,n 
^KK' 


=0{h^ 




MK G % 


VK' G Ni(K) 


Rl> n 


=0{h r,i 


+ 5tf) 


MK G % 


Me e£i(K)u£ w (K) 


F i,n 
K 


=0{K^ 


+ st?) 


MK G % 




-jyi,n 


=0{h r,a 


+ 6tf) 


Me G £ 




U i,n 


=0(K* 


+ 5tf) 


Ve G f 





We first prove an estimate for (I23p - (l24ft not taking into account the interface conditions: 



Lemma 4.3. Let {p l } u l {p 1 )) satisfy H23I) -(|24f) and suppose that assumptions 1 3. 1\ and \4.1 
are satisfied. 
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Then : 

1 2 2 Ni-1 , i,n+l _ i,n+l\2 

^EE(^) 2 ™-m + ±EE«.E E ^(m'l) 

2 JVi-1 ^> n+1 _ ^'"+^2 

+EE«-E E J) -~w 

i=l n=0 Ke% e£S t (K)U£ iD (K) y ' 

2 ^2 

< 2 ( n *'^)i2(0,T;L2(r)) + 2 E E (PK) 2meaS ( K ) 

i=i i=i A'eTi 

2 ATj-1 

i=l n=0 XgTj 
1 2 JVi-1 

+ 5 E E ^ E E &™<P% ¥1 - p£ +1 )meas{[KK>]) 

i=l n=0 AeTi K'^Ni{K) 
2 JVi-1 

" E E ^ E E Ri' n+1 pf +1 meas(s) 

i=l n=0 KeTi £G£ i (A)U£' lD (ii') 
2 JVj-1 

"EE^E E ^ n+1 (Pr +1 -^ n+1 )mea S ( £ ) (25) 

i=l n=0 eef^/Ou^i?) 

Proof. In each subdomain Oj, we multiply (|23p by b~tip l £ +1 and we sum over cells and 
time step n and make use of the following formula: 

E, i,n+l i,n\ i,n+l A X ■> \ / i,n+l\2 / i,n\2 , / i,n+l i,n\2 
(Pk ~Pk)Pk =2 L [(^ ) -{pkr + (pk -Pk) 

n=0 n=0 



1 

~2 
1 

^2 



n=0 



^-L d(x K ,x K ,y K u J7 2 ^ d(x K ,x K i) u J; 

Using the fact that Rrk' + Rk'k = (Assumption 14. lj) , we have : 

E E ^^'^ rimeas ([ i ^^ / ]) = E E R Tk'(pT' ~P l K )meas([KK']) 

KeTi K'dMi{K) K£Ti K'eAfi(K) 



16 



Then, we get: 

-> i Ni — l i i,n+l i,n+l\2 

\ E " (^°) 2 ) -as(if) + i J> £ £ df ^as^if']) 

Ni-1 Ni—l 
- E <^ E E nr +1 Px n+1 meas(e) < E <^ E F^+V^'meas^) 

n=0 A'eTi eefiC^UfinCK) n=0 Ke^ 

+ E E «x +1 ^ n+1 mea S (^)+i j>*E E ^(P^ +1 ^ +1 )nieas([^]) 
n=0 i^eTi n=0 iCe^ K'eN~i(K) 

Ni-1 

- E ^ E E i?r +1 ^ +1 mea S ( £ ) 

n=0 KsTi ee£i(K)U£ iD (K) 

We sum over the subdomains and make use of the following formula for any e £ £ j U 
and if = if (s), 

-<-V> = - nTvT - uTipt - pT) 

(i,n i,n \ 



i i,n t,n\2 
i,n i,n , v^K ~ -P £ j / i,n J.,n\jji,n 



Then, ([25]) follows from pl' n = for all e 6 £j£>. □ 

We focus now on the interface terms on T that appear in the first term of the right 
hand side of (|25p . We first consider the interface matching conditions ( ISi' ). 



Lemma 4.4. Let (p,u) satisfy ([25]) . ([2"4]). ([IS?]). TTien, we have: 

2 

E ( U ^).L 2 (0,T;.L 2 (r)) = ^ m '^)i 2 (0,T;L2(r)) + ^)L 2 (0,T;L 2 (r)) ( 2 6) 

1=1 
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Proof. From the interface conditions ( ISi' ), we have: 

N m -l N e -l 

St m u™' n+1 p™' n+1 meas(e) + 5t e ^ uf n+1 p e E ' n+1 meas(e) 

n=0 e££ n=0 ee£ 

= (w m ,p m ) L 2( 0)T;L 2( r )) + (U e ,p e ) L 2( 0jT;L 2( r )) 

= {Q m {~U e ) +^,p m ) L 2 (0iT;L 2 (r)) + (u e ,Q e ( P m ) +r S ) L 2 {0>T . L 2 {T)) 

= (-QmK),p m ) i2(0)T;L 2 {r)) + {U e , Qe(p m )> L 2 (0iT;L 2 {r)) 

v „ ' 

=0 (Lemma O 

+ (^£)P m )L 2 (0,T;L 2 (r)) + ('" e '^)L 2 (0,T;L 2 (r)) 



We consider now the interface scheme (IS^ ) 



Lemma 4.5. Let (p,u) satisfy (J23D, fl2U), PS 2 'D - Then, we have 



□ 



L 2 (0,T;L 2 (r)) 



< (p m ,^) L 2 (0)T;L 2 (r)) + K,^ " d £ U§ - d?(U £ - ^ m )>L 2 (0,T; 



L 2 (0,T;L 2 (r)) 



i=l 



(27) 



where d £ is the piecewise constant function on T such that d £ (x) = d\ for all x £ e 6 f . 
Proof. Using the interface conditions ( |IS2'D , we have: 

X^=l (^' ui )L 2 (0T-L 2 (r)) = (~Qm(u e ) +U£,P m )L 2 (0,T;L 2 (r)) 

+ (n e , Q e ( P m - d?u? + d^l^) + V £ - dfu% - d £ U£) L2(0jT . L2[T)) 
Since we use I? projection, we have 

Yh=1 {p % ' U% ) L 2 {0,T;L 2 {T)) = ( _ ^ + ^ £ ' P™) L 2 (0,T;L 2 (T)) 

+ {u ,p - d £ u £ +d £ U £ +Vs- d £ u £ - d £ U £ ) L2 ^ T . L2(X y } 



Simplifying the relation and using again the first equation of ([IS2J, 

Ei=l {p 1 i u1 )l 2 {0,T;L 2 (T)) = (U£'P m )L 2 (0,T;L 2 (r)) 

+ « -d?Q m {-u e ) - d?U £ + d?U? + V £ - dfu% - d £ U§) L2{0T . L2{r)) 
Since Q m is a L 2 projection, we have 

Yh=\ {P % ' U% ) L 2 {0T-L 2 (V)) - ^'-P m )L 2 (0,T;L 2 (r)) 

+ (u*, -dfU £ + dgU? + V £ - d £ U§) L2i0 T , LHr)) 



□ 



18 



Theorem 4.6. Suppose assumptions \ 3. 11 \4-l\ and \4-£\ hold. Let (p,u) satisfy (j23|) . ([Ml) , 

If one of the two conditions is satisfied, 



i) (p,u) satisfies transmission conditions fllSj'D , 



2P 



ii) (p,u) satisfies transmission conditions (IS2M and assumption \3.6\ holds 
Then, we have the following estimate: 

E \P%%,su +2E 11/^11^(00 < ll^°lli 2 ( 0i) + Oihf + 0(6t 

i=l i=l i=l 

where r\ = 2nuxi(j]j)j=i ! ... J 7 et 7 = 2min(7j)j=i ) ... ! 7. 

Proof. The proof consists in estimating the terms in formula (|25p of lemma 14.31 We often 
use the relation 

\ab\ ^ -a 2 + ^-b 2 V(a, b) € R 2 VC £ R* 
2 20 

with various constants C{ which are independent of the parameters of the mesh size. For 
the terms that are classical in the finite volume theory, we simply write the estimate. 
We begin with the estimate of the classical term 

2 Ni-l r 2 Ni-1 _ 

E Z> E ^ n+ V/ +1 meas(K) < f EE^E (pT +1 ) 2 ™™(K) + 0(h»* + 6$*) 

i=l n =0 K&Ti i=l n=0 K&Ti 

We also have the term 

2 Ni-l ^ 2 Ni-1 

E E E ^V^measW ^yEE^E (^ +1 ) 2 meas(^)+0(/ i ^+^) 

i=l n =0 ifeT; i=l n=0 ifeT; 

and the term 

2 JVj-1 „ 2 iVj-1 

-EE^E E K' n+1 PT +1 ™^) <y EE^E E (^' n+1 ) 2meas ( £ ) 

i=l n=0 KeT l£ eS t (K) i=l n=0 K&T ie e£i{K) 



+ 0(/i 2 " 4 + 5t 2li ) 
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We consider now 

2 Ni-l _ 

" £ Z> E £ ^' n+1 " Px n+1 )^ n+1 mea S (.) 

i=l n=0 K&% £ ^Si{K)Vj£ iD (K) 

2 JVi-1 i, n +l _ 

= - E E *** E E ^ lA( Ur +1 Vdt^T)me aS (e) 

i=l n=0 Ke% s££i{K)U£ iD (K) V d \ X K,ys) 

4ee^e E (ft ^ ' "-w 

1 2 

+ 2^E E ^ E E (^ +1 )M^,y £ )mea S (e) 

1 i=l n =0 A'eTi ee^(K)u£ iD (X) 

4ee*e e (Pe ^ - ' e > 

i=i n=o xer, e£«i(Jf)u£iD(]f) v ' ye; 

-. 2 iVi-1 

+ ^-O^ + ^ 71 )£ E ^ E E d(^,y e )meas( £ ) 

1 i=i n=0 AeTi eefi^UfioCA") 

„ 2 — 1 / j,n+l _ i,n+l\2 

4 E E * E E (P£ d , } + + «?* ) 

i=l n=0 XgT, £ g£ lWU £ lD (A) 

where we have used the following formula (see [7]) 

d{xKi X K') meas([K, K'}+ d{xx, He) meas(e) < d meas{Q,i) 

K£% K'£j\fi(K) e€£ t (K)U£ iD (K) 

In a classical way, we get 

, 2 Ni-l 

2 E £ E E - ^ +1 )meas([K^]) 

i=l n =0 ifeT; K'eAfi(K) 

n 2 Ni-l , i, n +l _ i,n+l\2 

4ee«.e e \: :)' -<w)^(>»» t ») 

1=1 n=0 AeT 4 K'eAfi(K) v ' ' 

and 

2 iVi-1 

E E E E - ^ n+1 )me as ( £ ) 

i=l n=0 XeT ieG £- i; (ii-)u£: i0 (if) 

„, 2 iVj— 1 / j,n+l _ j,n+l\2 

4 E E E E (P£ d , } + °(^ 4 + 
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We now focus on the interface terms in (|25|) starting with transmission scheme ( ISi'D . By 
Q26J) of lemma S3]: 

2 

E (u % ,P % ) r,2(n,r; r,2(r)) = (P m >^)L 2 (o,T ; L 2 (r)) + ^)i 2 (o,T ; L 2 (r)) 
i=i 

iV e — 1 e,n+l _ e,n+l 
+ y 5t e V^- ^_7^ +1 meas(e) 

JV m — 1 „ N e —1 , e,n+l _ e,n+l\2 

<f E ^E(^ +1 ) 2 meas( £ ) + ^ £ df ^ 

n=0 eG£ n=0 £ e£ V K,Ue) 

+ 0(/i) min( - ,?1+r ' 6 ' 2??6 ' 2,?7 - 1 + O (5* 2 ) min +T6 ' 276 ' 2 ^ 7 



The analysis of the interface term in (|25p with the transmission scheme (IS2M is more 
involved. By (|27p of lemma l3~o"l we have the following additional term: 

iVe-l 

- (u e ,dftf| + - Un)mo,T;L*(D) = ~ E ^ E ^(^ e ' n ) 2 meas(e) 

n=0 e€£e(Jf) 
v v ' 

=0(h 2 "i +5tl 11 ) 

Ne-1 p e > n _ p e > n iVe-l 

- E ^ E < ' d / e(£ V £ e -"meas( £ ) - ^ £ C^ ,n (W £ " - ^)meas( £ ) 

n=0 e££ e (K) 6 n=0 e££ e (K) 



N e — 1 e ' n e ' n 

-E^ E ~ I Kei£) m - t/r' n )meas( £ ) 

n=0 £ e£ e (K) £ 



/ e,n_ e,n »2 

By assumption l3.6l the last term is under control since ^^io 1 ^ e Yl £ e£ e (K) ^d? m eas(g) ^ 
aTmeas(0 m ). 

Summing up all these estimates, we have: 



^EE(^) 2 m eas(K )+ (l-^)EE^E E , {XK ll) meas{[KK ' ] 

fl_^_^_^_^_^W y y ^ ^ ~ P * } meas( g ) 

V 2 2 2 2 2 J ^ ^ ^ ^ d(x K ,y £ ) w 

V / i=l n =0 KeT, ee£ l (K)u£, D (K) V A ' ye; 
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2 „ „ 2 JVj-1 



i=l KeTi i=l n=0 ifGT; 

+ E E ^ E E (pr +1 ) 2 me as(£ ) 



i=l n=0 KGTi e&Si(K) 

n -\-i h . ) , / i , ^ 



+ 0(/i) min ^ 2 ^ i=1 --- 7,7?1+??8 ) + 0(<W2 s,min( - < - 27j ' >J=1 "' 7 ' 7l+76 ' > 



where C{ and C' 2 are zero for the scheme ( ISi' ) In other words, we have: 



1 X — ' II i,N l \\2 . • / 1 C3 n Cl C( C2 C4 ^6^^, i,2 

2 E lb 11^(00 +mm ( 2 ~ T' 1 2 2 2~"TJ^ Ip ll - r -^ 

i=i v 7 i=i 

^ 1 11 i,o 1 1 2 . C2 + C5 11 i\\2 C4 + C7 ^ .. j ||2 

^ 2 Z> l|p Hi 2 (0 l ) H 2 ^ " P l| i 2 (o^;i 2 (0 ! )) H 2 ^ " P lli2 (o,r ; L 2 (r)) 

i=l i=l i=l 

+ 0(/t) min ^ 2 ^ 3 ' =1 --- 7 '' n+ ' 78 - ) + 0(<W2) min( - < ' 27j ' >J=1 "' 7 ' 7l+76 ' > 

We notice that min(2r/i, 2r/e, ?7i + %) = 2min(r/i,ry6)- Moreover, using a result in [7] 
(discrete Poincare inequality, lemma 3.1), we have: 

HP i ||L»(o I T;L»(0 i )) ^diam(0 l )b 4 |i,r I ,5t 1 ViG {1,2} (28) 
We also use the discrete trace estimate proved in [8] : 

IIP'lllw^W) < C ( ^) (llfi^(o,T;L 2 (0O) + IPilr il5ti ) Vt e {1,2} (29) 

Taking small enough constants (independently of h) suffices to end the proof of Theo- 
rem WM n 



5 Well posedness 

Theorem 5.1. We assume that assumption \3. 1\ holds. 



Then, the problem defined by (|16l) - (fT5j) and interface scheme either fllSiD or QIS2D is well- 
posed. 

Proof. In both cases, we have a square linear system. It is thus sufficient to prove that the 
only solution with a zero right hand side and initial condition is zero, i.e. 

U* e < n = = #£ K , = Rl> n = F^ n =p l £ = V £ =U £ = ViG {1,2} VnG{0...A^} 
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Then, by Lemma l4.3l and Lemma l4.4l (resp. 14. 5p for transmission scheme pSiD (resp. QIS2D 
we have 



, 2 1 2 Ni—l 



(p^ +1 -^ n+1 ) 2 



i=i /reTi 



i=l n=0 Ke^ K'eMi(K) \ A ; i\ 



meas( [i^iT']) 



2 JVi-l 

+EE«-E E 

i=l n =0 XgT, ee£i(A")U£ iD (i ; £') 



(pW j,n+l )2 

^ meas(e) < 

d{XK,ye) 



that is for alH = 1, 2, 1 < n < iVj, p l ' n is a constant. By equation (|16p in every subdomain, 
we then have that the value of the constant is independent of n. Since the initial condition 
is zero, the constant is actually zero. □ 



6 Error estimate 



Let p , p 2 be the solution to the continuous problem (I13al) - (|13d[) . We define the interpola- 
tion on the mesh T\ U T2 at time t n by: 





= p l (x K ,t n ) 




1 f*n+l/2 


vT 


6ti k-i/2 


vT 


= Vee£ iD 




1 f t n + l/2 




St, Li 

1 Jt n -l/2 



meas e 



p i Ve G Si 



^ Vee£iU£ iD 



> VnG{0,...,iVi} ViG {1,2} 



meas(e) J £ dn 

We have to estimate the error terms e l K = p l K — p l K , e\ = p\ — p\ and q\ = u\ — u\ . 
Theorem 6.1. We suppose that the solution has the following regularity: 

P eC\O,T;C 2 (0)) (30) 

and that the numerical right hand side is such that: 
1 rti 



TT I ^ / -f)= 0{diam{K) + SU) Vn G {0, 

dti J tl meas(K) J K V " J 

1/2 



,Ni} Vi€{l,2} 
(31) 



We assume that Assumptions \3. 11 3.7 hold and that transmission scheme (ISi) is used. 
Then, we have the following estimate: 

2 2 

E + E 11^*11^(00 = 0(h + St) 



1=1 



1=1 
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where 5t = max(5ti, ^2). 

The same estimate holds for transmission scheme QIS2D if in addition Assumption ^. 61 holds. 



Proof. It is easy to check that the errors ex, e £ et q £ satisfy (l23]) - (l24j) with error terms 
defined by 

„i,n = J_ f tll+1 / 2 1 f ( A,n _ f 
K StJ tii/2 me aS (K)J K V K 1 



-t ft} 1 f r> i ~i,n ~i,n—l 

1 fn+1/2 1 r dp Px ~Pk 



5t{ J t i i 2 meas(if) J K dt 5t 

R i,n J_ /''™+ 1 /2 1 f d]f__ p l (x K ,,t n ) -p l (x K ,t n ) 

KK ' 5tiJ e meas([KK'}) J< KK n dni d{x K ,x K >) 

n — 1/2 L J 

dU J t i meas(e) J e dm 

n — 1/2 



d>( x Ki(e),Ue) 
Ke,T = Qm{-U e £yT ) - u™ T 

For the scheme (|TS7]): P £) r = Qe{p? tT ) ~ Pt,T 



. For the scheme flE^: 7> £ , T = Q e ^ m(e) J - p^ e(e) - (C + t£)££ T 

The derivation of the formula for V Et T and U E) T are made explicit in the sequel when these 
terms are estimated. Let us remark that we have by construction e l £ = for all K E 
i = 1,2. 

By assumption and by using Taylor expansion, it is classical to check that the error 
terms Rk, Rkk' and Fk satisfy Assumption 14.21 with rji = 7i = 1 for i = 2, 3, 5. As regards 
the term U e , we proceed as in pQ. From Assumption (|30p on the regularity of the solution, 
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we have: 



jji,n " i - / / ' 



d(x K j £ ),y £ ) Sti J t i meas(e) J £ drii 

^^C^^k^^P 1 -P^Ve^n) ^ f p i (y £ ,t n )-p i (x Ki{£) ,t n ) dp* \ 

d(x Kl (e),ye) \ d(x K .^,y £ ) drii £ ' n J 

( dp* t J_ 1 r dp*\ 

[drii e ' 5U J t i meas(e) J £ drii J 

\ n — 1/2 / 

0(diam(e)) 2 

=— + 0(d(xx i ( £ ),ye)) + 0(diam(e)) = O(h) by assumption 13.71 

d\XKi(e)i He) 

thus, r\\ = 1 and 71 > 1. Since, 

6ti J t * 2 meas(e) J £ ydrii 5U J? 2 meas(e) J £ drii J 
we have 774, 74 > 1. 

We now consider the non classical consistency error terms V £ et IA £ . For the transmission 



condition ( |ISiD we have: 

q?=u?-uf = Q m {-u^-u? 



We 

The error on the transmission condition on the interface reads: 

e e £ =p e £ -p e £ = Qe(p?)-p e £ 

= Q e (e?+p?)-p e £ = Q e (e?) + Qe{p?)-P e £ 

If (m, e) = (1, 2), we get : 

K = [Q2(pl, T )] n -fe' n 

5t2 fc^i 5tl Jt n ,k-i/2 meas ( e ) A ^2 7i n _ 1/2 meas(e) J £ 



U n,k = ^ Qi{ _~2 ) n,k_ {i l 



n, k 



t n 1 ,fe+l/2 



1 /-*«+i/2 1 r dp 2 1 /•*» 1 rap 1 , 



5t2 Jt n _ 1/2 m eas(e) j e 9n2 5ii Jt k n ~ Yl2 meas(e) j e dri\ 
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If (m, e) = (2, 1), we have 



V. 



n,k 



[QM,T)r k -Pi' n ' k 



1 



*n+l/2 



1 



8to h , meas(e) 

z Jt n-l/2 y ' 

21 n 



P 



St- ^ £ 



el u el 

K. 



k=l 



k + l/2 



^2 ~ <fti A*-V 2 meas(e) 7 e <9ni 



dp 1 



fc+l/2 



1 



1 



Jt k ~ 1/2 meas(e) 



*n + l/2 



St 2 Jt n _ 1/2 incas( -) 



9ng 







We consider now transmission scheme (IS2): IA £ is left unchanged and V £ now reads 

Q e ( p ™-d? U ?)-pt-d?ut 
Q e (e? + " ~ Cu7) -Pt- d ?u e £ 

Q e (e? - d?q? + dfUT) - d%Ut + Q e fe (£) ) - p e Ke{£) - (df + dt)ul 



If (m, e) = (1, 2), we have 



^ = [Q2(^ l(£) )] n -4: ( e)-(^+^ 



;2\~2,n 



<5t 1 f 

-^Y^P l {x Kl{e) ,t k n ) -p 2 (x K2{e) ,t n ) - (d\ + d 2 £ )— J _ 



1+1/2 



fc=l 

a: 



; meas(e) j £ <9ri2 

n— 1/2 



dp 2 



— y 



fe+1/2 



2— ^tr i/2 meas(e ;Je 
=0(/t + 5t) 
If (m, e) = (2, 1), V £ reads : 



^(^(e),^) ~ P 2 ( x K 2 (e),tn) - d{x Kl (e),X K2 (e)) 



dp 2 
dm 



n,k 



1 



2 

K 

,fe+l/2 



2i (Pjc a ( e ) 



n,k 



~l,n,k , ,1 . t2\~ 
Pjfi(e) - ( d e + d e) u . 



fe-1/2 meas e 



P 2 OftT 2 (e) ' *n ) ~ P 1 ( x ^i (e) » *n) ~ d ( x Xx (e) , ^ (e) ) 



dp 1 

9ni 



=0(/i + 5t) 
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For all cases, we have thus ?76> ??7) 76> 77 > 1- Then the error estimate follows from Theo- 
rem HjJJ □ 
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7 Numerical results 



In this part, we illustrate the method with a parabolic equation coming from a previous 
article of Ewing and Lazarov We consider the (IS2) interface conditions, i.e. equation 
(|l(jp in one dimension, which are more natural. We solve the following model problem : 



^-p( x ,t) - ^(x,t) = f(x,t) Vte[0.,0.1] VxG[0.,l.] 



p(x,t) 
p(x,0) 



0,Vx G dtt 




(32) 

(33) 
(34) 



The following function is used as an exact solution : 

p(x, t) = exp(20(t - t 2 ) - 37x 2 + 8x - 1) 
This function represents a bump with a maximum value near the position x = 0.15. 




Figure 4: Exact solution for t=0.1 

In the interval [0.5,1.], the function is close to 0. In this interval, the function changes 
negligibly in time. In contrast, the function changes rapidly in time in the interval [0., 0.5] 
and simulates a local behavior. 

We use two different time step sizes : 

• a fine time step 5t\ = 0.002 in the subdomain VL\ = [0.,0.25], discretized with a fine 
grid 5xi = 0.01 

• a coarse time step 5t2 = 0.02 in the subdomain VI2 = [0.25,1.], discretized with a 
coarse grid 6x2 = 0.05 

The interface is placed at x = 0.25 . This is a worst case since the domain with local 
refinement only partially covers the interval [0.,0.25] where the solution changes quickly. 
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We consider two cases. The first one (coarse master) is when the coarse domain enforces 
the Dirichlet condition, see equation (|l(jp . The second one (fine master) is when the refined 
domain enforces the Dirichlet condition, see equation (jlip . We make a comparison with 
the algorithm given by Mlacnik and Heinemann [9j [10] . In the following pictures, we plot 
the evolution of the errors in space and the time evolution of the L 2 norm of the error 
between the exact solution, the two local time step methods and the solution with the fine 
or coarse time step on the whole domain. 

At each coarse time step, we solve the set of discretized equations using the iterative 
algorithm explained in section [2T21 with the stopping criterium e = 10~ 5 . In the following 
figures, we plot the error between these two solutions and the exact solution. For com- 
pleteness, we also plot the error for a computation with either the coarse or the fine time 
step on the whole domain. The number of iterations needed to reach the convergence is 
quite small; it is about 6 for the fine master method and about 8 for the coarse master 
method. 

We notice that for both cases, the error is significantly smaller than the one of the 
coarse time step. Morever, for the fine master method, see equation tjlip . the error is close 
to the fine time step error in the refined zone, see figure [5j 

As explained in section [2T2l it is not necessary to iterate until convergence the algorithm 
to obtain a conservative method. In figure (U we plot the error after only one iteration of 
the corrector stage. As expected, the errors are larger than with the converged solutions. 
Let us recall, that the method proposed by Mlacnik corresponds to the fine master curve 
in figure El 
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Figure 5: After convergence, on left : error in space with the exact solution at time t = 
on right : time evolution of the 1? norm of the error. 



/ \ 


'coarse dt' — i — 

fine masler * 
coarse master □ 


" / \ 

/ \ 




1^' 1 : : C" 


.....-B B--..43 *"""5-. ^ 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.02 0.03 0.04 0.05 0.06 0.07 0.06 0.09 0.1 0.11 

Figure 6: After one iteration, on left : error in space with the exact solution at time t = 
on right : time evolution of the 1? norm of the error. 
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8 Conclusion 



We have proposed a local time step strategy for solving problems on grids with different 
time steps in different regions. We have analyzed two schemes: QISiD an d PS2D • I n PSiD, 



the coupling involves additional interface unknowns. Scheme ( IS2 ) is written only in terms 
of "classical" finite volume unknowns. Both schemes are conservative, of order one in 



space and time. The assumptions are more restrictive for ( IIS2D than for JISi [ ) . We have 
presented an iterative solution method for solving the composite grid system. Its main 
feature is that at every stage, conservativity is ensured. Numerical tests on a toy problem 
confirm the capabilities of the method. The scheme is being implemented in a multiphase 
three dimensional simulation code. 
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